#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <float.h>
#include <limits.h>
#include <assert.h>
#include <mpi.h>

#include "halofinder.h"

int read_GADGETsinglefile_Page(long PageNum, long *NumBuffParts, float **buffPx, float **buffPy, float **buffPz, float **buffVx, float **buffVy, float **buffVz)
{
  static FILE *fp;
  long k;
  char basefname[MAX_FILENAME];
  char fname[MAX_FILENAME];
  float *partChunk,sqrta;
  long TotNumPageParts;
  long NumPartsInChunk,BaseNumPartsInChunk,NumChunks;
  long seekval;
  fpos_t pos_p,pos_v;

  struct io_header_1 {
    unsigned int npart[6];
    double   mass[6];
    double   time;
    double   redshift;
    int      flag_sfr;
    int      flag_feedback;
    unsigned int      npartTotal[6];
    int      flag_cooling;
    int      num_files;
    double   BoxSize;
    double   Omega0;
    double   OmegaLambda;
    double   HubbleParam;
    char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  // fills to 256 Bytes
  } header1;
  
  int dummy;
  
#define SKIP fread(&dummy,sizeof(dummy),(size_t) 1,fp);
  
  fp = fopen(haloFinderData.SnapshotFileList,"r");
  assert(fp != NULL);
  fgets(basefname,MAX_FILENAME,fp);
  fclose(fp);
  
  fp = fopen(basefname,"r");
  if(fp == NULL)
    {
      sprintf(fname,"%s.%d",basefname,(int) PageNum);
      fp = fopen(fname,"r");
      assert(fp != NULL);
    }
  
  SKIP;
  fread(&header1,sizeof(header1),(size_t) 1,fp);
  SKIP;
  
  TotNumPageParts = 0;
  for(k=0;k<5;k++)
    TotNumPageParts += (long) (header1.npart[k]);
  
  sqrta = (float) (sqrt(header1.time));

  BaseNumPartsInChunk = (long) (IO_CHUNK_SIZE_MB*1024.0*1024.0/((double) (sizeof(float)))/3.0); //IO_CHUNK_SIZE_MB MB read chunks for parts
  partChunk = (float*)malloc(sizeof(float)*3*BaseNumPartsInChunk);
  assert(partChunk != NULL);
  NumChunks = TotNumPageParts/BaseNumPartsInChunk;
  if(NumChunks*BaseNumPartsInChunk != TotNumPageParts)
    NumChunks += 1;
  
#ifdef DEBUG
#if DEBUG_LEVEL > 1      
  fprintf(stderr,"%d: NumPart = %ld, NumPartsInChunk = %ld, NumChunks = %ld \n",ThisTask,TotNumPageParts,
	  BaseNumPartsInChunk,NumChunks);
  fprintf(stderr,"%d: header1.npart[1] = %u \n",ThisTask,header1.npart[1]);
#endif
#endif
  
  //get position of velocity and particle records
  SKIP;
  fgetpos(fp,&pos_p);
  fseek(fp,(long) dummy,SEEK_CUR);
  SKIP;
      
  SKIP;
  fgetpos(fp,&pos_v);
  
  //sets how many bytes to seek to proper chunk
  seekval = PageNum*sizeof(float)*3*BaseNumPartsInChunk;
  
  //gets actual # of parts in chunk
  if(PageNum == NumChunks-1)
    NumPartsInChunk =  TotNumPageParts - PageNum*BaseNumPartsInChunk;
  else
    NumPartsInChunk = BaseNumPartsInChunk;
  
  //alloc mem
  *NumBuffParts = NumPartsInChunk;
  *buffPx = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffPx) != NULL);
  *buffPy = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffPy) != NULL);
  *buffPz = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffPz) != NULL);
  *buffVx = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffVx) != NULL);
  *buffVy = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffVy) != NULL);
  *buffVz = (float*)malloc(sizeof(float)*NumPartsInChunk);
  assert((*buffVz) != NULL);
  
  //read positions
  fsetpos(fp,&pos_p);
  fseek(fp,seekval,SEEK_CUR);
  fread(partChunk,(size_t) NumPartsInChunk,3*sizeof(float),fp);
  for(k=0;k<NumPartsInChunk;++k)
    {
      (*buffPx)[k] = partChunk[k*3 + 0];
      (*buffPy)[k] = partChunk[k*3 + 1];
      (*buffPz)[k] = partChunk[k*3 + 2];
    }
  
  //read velocities
  fsetpos(fp,&pos_v);
  fseek(fp,seekval,SEEK_CUR);
  fread(partChunk,(size_t) NumPartsInChunk,3*sizeof(float),fp);
  for(k=0;k<NumPartsInChunk;++k)
    {
      (*buffVx)[k] = (float) (partChunk[k*3 + 0]/sqrta);
      (*buffVy)[k] = (float) (partChunk[k*3 + 1]/sqrta);
      (*buffVz)[k] = (float) (partChunk[k*3 + 2]/sqrta);
    }
  
  //clean up and return
  fclose(fp);
  free(partChunk);
  
  return 1;
      
#undef SKIP  
}

long read_GADGETsinglefile_NumPages(void)
{
  FILE *fp;
  char basefname[MAX_FILENAME];
  char fname[MAX_FILENAME];
  long TotNumPageParts,k;
  long BaseNumPartsInChunk,NumChunks;
  
  struct io_header_1
   {
     unsigned int      npart[6];
     double   mass[6];
     double   time;
     double   redshift;
     int      flag_sfr;
     int      flag_feedback;
     unsigned int      npartTotal[6];
     int      flag_cooling;
     int      num_files;
     double   BoxSize;
     double   Omega0;
     double   OmegaLambda;
     double   HubbleParam;
     char     fill[256- 6*4- 6*8- 2*8- 2*4- 6*4- 2*4 - 4*8];  /* fills to 256 Bytes */
   } header1;
  int dummy;
  
#define SKIP fread(&dummy,sizeof(dummy),(size_t) 1,fp);
   
   fp = fopen(haloFinderData.SnapshotFileList,"r");
   assert(fp != NULL);
   fgets(basefname,MAX_FILENAME,fp);
   fclose(fp);
   
#ifdef DEBUG
#if DEBUG_LEVEL > 1
   fprintf(stderr,"%d: reading num pages from '%s'\n",ThisTask,basefname);
#endif
#endif
   fp = fopen(basefname,"r");
   if(fp == NULL)
     {
       sprintf(fname,"%s.%d",basefname,0);
       fp = fopen(fname,"r");
       assert(fp != NULL);
     }
   
   SKIP;
   fread(&header1,sizeof(header1),(size_t) 1,fp);
   SKIP;
   
   fclose(fp);
   
   TotNumPageParts = 0;
   for(k=0;k<5;k++)
     TotNumPageParts += (long) (header1.npart[k]);
   
   BaseNumPartsInChunk = (long) (IO_CHUNK_SIZE_MB*1024.0*1024.0/((double) (sizeof(float)))/3.0); //IO_CHUNK_SIZE_MB MB read chunks for parts
   NumChunks = TotNumPageParts/BaseNumPartsInChunk;
   if(NumChunks*BaseNumPartsInChunk != TotNumPageParts)
     NumChunks += 1;
      
#ifdef DEBUG
#if DEBUG_LEVEL > 1      
   fprintf(stderr,"%d: NumPart = %ld, NumPartsInChunk = %ld, NumChunks = %ld \n",ThisTask,TotNumPageParts,
	   BaseNumPartsInChunk,NumChunks);
   fprintf(stderr,"%d: header1.npart[1] = %u \n",ThisTask,header1.npart[1]);
#endif
#endif
   
#undef SKIP  
   
   return NumChunks;
}

